Phase mixing of shear Alfven waves as a new mechanism for electron acceleration in 

collisionless, kinetic plasmas 
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Particle-In-Cell (kinetic) simulations of shear Alfven wave interaction with one dimensional, across 
the uniform magnetic field, density inhomogeneity (phase mixing) in collisionless plasma were per- 
formed for the first time. As a result a new electron acceleration mechanism is discovered. Progres- 
sive distortion of Alfven wave front, due to the differences in local Alfven speed, generates nearly 
parallel to the magnetic field electrostatic fields, which accelerate electrons via Landau damping. 
Surprisingly, amplitude decay law in the inhomogeneous regions, in the kinetic regime, is the same 
as in the MHD approximation described by Heyvaerts and Priest (1983). 

PACS numbers: 52.65.Rr; 52.35.Bj; 52.35.Mw; 96.60.Ly 



Interaction of shear Alfven waves (AWs) with plasma 
inhomogeneities is a topic of considerable importance 
both in astrophysical and laboratory plasmas. This is 
due to the fact both AWs and inhomogeneities often co- 
exist in many of these physical systems, shear AWs are 
believed to be good candidates for plasma heating, en- 
ergy and momentum transport. On the one hand, in 
many physical situations AWs are easily excitable and 
they are present in a number of astrophysical systems. 
On the other hand, these waves dissipate on shear vis- 
cosity as opposed to compressive fast and slow magne- 
tosonic waves which dissipate on bulk viscosity. In as- 
trophysical plasmas shear viscosity is extremely small as 
compared to bulk one. Hence, AWs are notoriously dif- 
ficult to dissipate. One of the possibilities to improve 
AW dissipation is to introduce progressively collapsing 
spatial scales, 81 — > 0, into the system (recall that the 
classical viscous and ohmic dissipation is oc 5l~ 2 ). Hey- 
vaerts and Priest have proposed (in astrophysical con- 
text) one such mechanism called AW phase mixing [l|. It 
occurs when a linearly polarised shear AW propagates in 
the plasma with one dimensional, transverse to the uni- 
form magnetic field density inhomogeneity. In such sit- 
uation initially plane AW front is progressively distorted 
because of different Alfven speeds across the field. This 
creates progressively strong gradients across the field (ef- 
fectively in the inhomogeneous regions transverse scale 
collapses), and thus in the case of finite resistivity, dis- 
sipation is greatly enhanced. Thus, it is believed that 
phase mixing can provide substantial plasma heating. A 
significant amount of work has been done in the context 
of heating open magnetic structures in the solar corona 
B H H H la 11 II HH ES El ■ All phase mixing studies 
so far have been performed in the MHD approximation, 
however, since the transverse scales in the AW collapse 
progressively to zero, MHD approximation is inevitably 
violated, first, when the transverse scale approaches ion 
gyro-radius n and then electron gyro-radius r e . Thus, 
we proposed to study phase mixing effect in the kinetic 



regime, i.e. we go beyond MHD approximation. As a 
result we discovered new mechanism of electron acceler- 
ation due to wave-particle interactions which has impor- 
tant implications for various space and laboratory plas- 
mas, e.g. coronal heating problem and acceleration of 
solar wind. 

We used 2D3V, the fully relativistic, electromag- 
netic, particle-in-cell (PIC) code with MPI parallelisa- 
tion, modified from 3D3V TRISTAN code 12]. The sys- 
tem size is L x = 5000A and L v = 200A, where A(= 1.0) 
is the grid size. The periodic boundary conditions for 
x- and y-directions are imposed on particles and fields. 
There are about 478 million electrons and ions in the 
simulation. The average number of particles per cell is 
100 in low density regions (see below). Thermal velocity 
of electrons is v t h. e = 0.1c and for ions is Vth,i — 0.025c. 
The ion to electron mass ratio is mi/m e — 16. The time 
step is uj pe At — 0.05. Here uj pe is the electron plasma 



frequency. The Debye length is v t h, e /u 



pe 



1.0. The 

electron skin depth is c/oj pe = 10 A, while the ion skin 
depth is c/ujpi — 40 A. Here cu p i is the ion plasma fre- 
quency. The electron Larmor radius is Vth.e/^ce = 1.0A, 
while the same for ions is Vth,i/^d = 4. OA. The external 
uniform magnetic field, Bq, is in the x-direction and the 
initial electric field is zero. The ratio of electron cyclotron 
frequency to electron plasma frequency is Lu ce /uj pe = 1.0, 
while the same for ions is uj C i/uj p i — 0.25. The latter 
ratio is essentially Va/c - the Alfven speed. Plasma 
[3 = 2(u} pe /uj ce ) 2 (vth,e/c) 2 — 0.02. Here all plasma pa- 
rameters are quoted far away from the density inhomo- 
geneity region. The dimensionless (normalised to some 
reference constant value of no) ion and electron density 
inhomogeneity is described by 



m(y) = n e (y) = 1+3 cxp - {(y - 100A)/(50A)) 



= F{y)- 

This means that in the central region (across y-direction) 
density is smoothly enhanced by a factor of 4, and 
there are strong density gradients of width of about 50A 



2 



around the points y = 51. 5A and y = 148. 5A. The 
background temperature of ions and electrons, and their 
thermal velocities are varied accordingly 

T l (y)/To = T e ( y )/To = F( y y\ (2) 



Vth,i/Vi0 = Vth,e/VeO = F(y) 1/2 , (3) 

such that the thermal pressure remains constant. Since 
the background magnetic field along x-coordinate is also 
constant, the total pressure remains constant too. Then 
we impose current of the following form 

3 t E y = -J„ sin(w«,t) (1 - exp [~(t/t ) 2 ]) , (4) 



d t E z = -J cos{cj d t) (1 - exp [-{t/t ) 2 ]) . (5) 

Here ua is the driving frequency which was fixed at u>d — 
0.3w c i, which ensures that no significant ion-cyclotron 
damping is present. Also, dt denotes time derivative, to 
is the onset time of the driver, which was fixed at 50/w pe 
or 3.125/w c j. This means that the driver onset time is 
about 3 ion-cyclotron periods. Imposing such current on 
the system results in the generation of left circularly po- 
larised shear AW, which is driven at the left boundary 
of simulation box and has a width of 1A. Initial ampli- 
tude of the current is such that relative AW amplitude is 
about 5 % of the background in the low density regions, 
thus the simulation is weakly non-linear. 

Because of no initial (perpendicular to the external 
magnetic field) velocity excitation was imposed in addi- 
tion to the above specified currents (cf. the excited 
(driven) at the left boundary circularly polarised AW is 
split into two circularly polarised AWs that travel in op- 
posite directions. The dynamics of these waves is shown 
in Fig.l where we show three snapshots of the evolution. 
Typical simulation, till the last snapshot shown in the fig- 
ure, takes about 8 days on the parallel 32 dual 2.4 GHz 
Xeon processors. It can be seen from the figure that be- 
cause of the periodic boundary conditions circularly po- 
larised AW that was travelling to the left has reappeared 
from the right side of the simulation box (t = 15.62/o; C j). 
Then the dynamics of the AW (B z ,E y ) progresses in 
a similar manner as in MHD, i.e. it phase mixes 0. 
In other words middle region (in y-coordinate) travels 
slower because of the density enhancement (note that 
Va (y) oc 1/ \/ni(y)) . This obviously causes distortion of 
initially plane wave front and the creation of strong gra- 
dients in the regions around y — 50 and 150. In MHD ap- 
proximation, in the case of finite resistivity 77, in these re- 
gions AW is strongly dissipated due to enhanced viscous 
and ohmic dissipation. This effectively means that the 
outer and inner parts of the travelling AW are detached 
from each other and propagate independently. This is 
why the effect is called phase mixing - after long time, in 
the case of developed phase mixing, phases in the wave 



front become effectively uncorrelated. A priori it was not 
clear what to expect from our PIC simulation because it 
was performed for the first time. The code is collision- 
less and there are no sources of dissipation in it (apart 
from possibility of wave-particle interactions). It is evi- 
dent from Fig.l that at later stages (t = 54.69/w C i) AW 
front is strongly damped in the strong density gradient 
regions. This immediately rises a question of where AW 
energy went? The answer lies in Fig. 2, where we plot E x 
longitudinal electrostatic field, and electron phase space 
(V x /c vs. x) for the different times (note, that in order 
to reduce figure size, only electrons with V x > 0.15c were 
plotted). In the regions around y = 50 and 150 for later 
times significant electrostatic field is generated. This is 
the consequence of stretching of AW front in those re- 
gions because of difference in local Alfven speed. In the 
right column of this figure we see that exactly in those 
regions where E x is generated, the electrons are acceler- 
ated in large numbers. Thus, we conclude that energy of 
the phase mixed AW goes into acceleration of electrons. 
In fact line plots of E x show that this electrostatic field 
is strongly damped, i.e. energy is channelled to electrons 
via Landau damping. 

The next piece of evidence comes from looking at the 
distribution function of electrons before and after the 
phase mixing took place. In Fig. 3 we plot distribution 
function of electrons at t = and t = 54.69/w c ;. Note 
that even at t = Q the distribution function does not look 
as purely Maxwellian because of the fact that tempera- 
ture varies across y— coordinate (to keep total pressure 
constant) and the graph is produced for the entire sim- 
ulation domain. There is also a substantial difference at 
t = 54.69/w C i to its original form because of the afore- 
mentioned electron acceleration. We see that the num- 
ber of electrons having velocities V x = ±(0.1 — 0.3)c is 
increased. Note that the acceleration of electrons takes 
place mostly along the external magnetic field (along x- 
coordinate). No electron acceleration occurs in V y or V z 
(not plotted here). 

The next step is to check whether the increase in elec- 
tron velocities comes from the resonant wave particle in- 
teractions. For this purpose in Fig. 4 we plot two snap- 
shots of Alfven wave B z (x,y — 148) components at in- 
stances t — 54.69/a; C i (solid line) and t — 46.87/w c , (dot- 
ted line). The distance between the two upper leftmost 
peaks (which is the distance travelled by the wave in time 
between the snapshots) is about 8L = 150A = 15(c/<x> pe )- 
Time difference between the snapshots is 5t = 7.82/u> C i. 
Thus, measured AW speed at the point of the strongest 
density gradient (y = 148) is Vf = SL/St = 0.12c. 
Now we can also work out the Alfven speed. In the 
homogeneous low density region the Alfven speed was 
set to be V^(oo) = 0.25c. From Eq.(l) it follows that 
for y — 148 density is increased by a factor of 2.37 
which means that the Alfven speed at this position is 
Va(148) = 0.25/V2T37c = 0.16c. The measured and cal- 
culated Alfven speeds in the inhomogeneous region do 
not coincide. This can be attributed to the fact that 
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in the inhomogeneous regions (where electron accelera- 
tion takes place) because of momentum conservation AW 
front is decelerated as it passes on energy and momen- 
tum to the electrons. However, this may be not the case 
if wave-particle interactions play the same role as dissipa- 
tion in the MHD: Then wave-particle interactions would 
result only in the decrease of the AW amplitude (dis- 
sipation) not in its deceleration. If we compare these 
values to Fig. 3, we deduce that these are the velocities 
> 0.12c above which electron numbers with higher veloc- 
ities are greatly increased. This deviation peaks at about 
0.25c which in fact corresponds to the Alfven speed in the 
lower density regions. This can be explained by the fact 
the electron acceleration takes place in wide regions (cf. 
Fig. 2) along and around y = 148 (and y — 51) - thus 
the spread in the accelerated velocities. 

In Fig. 4 we also plot a visual fit curve (dashed line) in 
order to quantify the amplitude decay law for the AW 
(at t = 54.69/a; C i) in the strongest density inhomogene- 
ity region. The fitted (dashed) cure is represented by 



0.056 exp 



(x/1250y 



There is an astonishing simi- 
larity of this fit with the MHD approximation results. 
Authors of Ref. found that for large times (devel- 
oped phase mixing), in the case of harmonic driver, the 

amplitude decay law is given by oc exp — ( ''gys' 4 ) x 3 

which is much faster than the usual resistivity dissipa- 
tion oc exp(— r]x). Here V' A is the derivative of the Alfven 
speed with respect to y-coordinate. The most intriguing 
fact is that even in the kinetic approximation the same 
oc exp(-Ac 3 ) law holds as in the MHD. In the MHD a 
finite resistivity and Alfven speed non-uniformity are re- 
sponsible for the enhanced dissipation via phase mixing 
mechanism. In our PIC simulations (kinetic phase mix- 
ing), however, we do not have dissipation and collisions 
(dissipation). Thus, in our case wave-particle interac- 
tions play the same role as resistivity r\ in the MHD phase 
mixing. It should be noted that no significant AW dissi- 
pation was found away from the density inhomogeneity 
regions. This has the same explanation as in the case 
of MHD - it is the regions of density of inhomogeneities 
(V A ^ 0) where the dissipation is greatly enhanced, while 
in the regions where V A = there is no substantial dissi- 
pation (apart from classical oc exp(— r/x) one). In MHD 
approximation aforementioned amplitude decay law is 
derived from the diffusion equation, to which MHD equa- 
tions reduce to for large times (developed phase mixing) . 
It seems that kinetic description leads to the same type of 
diffusion equation. It is unclear although, at this stage, 
what physical quantity would play role of resistivity r\ 
(from the MHD approximation) in the kinetic regime. 

It is worthwhile to mention that in MHD approxima- 
tion authors of Refs. [H, [Tl| showed that in the case of 
localised Alfven pulses, Heyvaerts and Priest's amplitude 
decay formula oc exp(—Ax 3 ) (which is true for harmonic 
AWs) is replaced by the power law B z oc x~ 3 l 2 . A nat- 
ural next step forward would be to check whether in the 
case of localised Alfven pulses the same power law holds 



in the kinetic regime. 

Finally we would like to mention that after this study 
was com plet e we became aware of study by Vasquez and 
Hollweg [l3|, who used hybrid code (electrons treated as 
neutralising fluid, while ion kinetics is retained) as op- 
posed to our (fully kinetic) PIC code, to simulate reso- 
nant absorption. They found that a planar (body) Alfven 
wave propagating at less than 90° to a background gradi- 
ent has field lines which lose wave energy to another set of 
field lines by cross- field transport. Further, Vasquez 0] 
found that when perpendicular scales of order 10 proton 
inertial lengths (lOc/ujpi) develop from wave refraction in 
the vicinity of the resonant field lines, a non-propagating 
density fluctuation begins to grow to large amplitudes. 
This saturates by exciting highly oblique, compressive, 
and low-frequency waves which dissipate and heat pro- 
tons. These processes lead to a faster development of 
small scales across the magnetic field, i.e. this is ulti- 
mately related to the phase mixing mechanism, studied 
here. Continuing this argument we would like make clear 
distinction between the effects of phase mixing and reso- 
nant absorption of the shear Alfven waves. Historically, 
there was some confusion in use of there terms. In fact, 
earlier studies which were performed in the context of 
laboratory plasmas, e.g. Ref. who proposed (quote) 
"the heating of collisionless plasma by utilising a spa- 
tial phase mixing by shear Alfven wave resonance and 
discussed potential applications to toroidal plasma" used 
term phase mixing to discuss the physical effect which 
in reality was the resonant absorption. In their later 
works, Hasegawa and Chen 0] and Hasegawa and Chen 
|17| . when they treated the same problem in the kinetic 
approximation (note that Hasegawa and Chen |15| used 
MHD approach) the authors avoided further use of term 
phase mixing and instead they talked about linear mode 
conversion of shear Alfven wave into kinetic Alfven wave 
near the resonance oj 2 — fcjj Va(x) (in their geometry the 
inhomogeneity was across the x-axis) . In the kinetic ap- 
proximation Hasegawa and Chen 16] and Hasegawa and 
Chen [ljfl found a number of interesting effects including: 

(i) They established that near the resonance initial 
shear Alfven wave is linearly converted into kinetic 
Alfven wave (which is the Alfven wave modified by the 
finite ion Larmor radius and electron inertia) that propa- 
gates almost parallel to the magnetic field into the higher 
density side. Inclusion of the finite ion Larmor radius 
and electron inertia removes the logarithmic singularity 
present in the MHD resonant absorption which exists be- 
cause in MHD shear Alfven wave cannot propagate across 
the magnetic field. Physically this means that finite ion 
Larmor radius prevents ting ions to magnetic field lines 
and thus allows propagation across the field. The elec- 
tron inertia eliminates singularity in the same fashion but 
on different, ri^/Te/Ti, spatial scale (here ion gyro- 
radius). 

(ii) They also found that in collionless regime (which 
is applicable to our case) dissipation of the (mode con- 
verted) kinetic Alfven wave is through Landau damping 
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of parallel electric fields both by electrons and ions. How- 
ever, in the low plasma (3 regime (which is also applica- 
ble to our case (3 = 2(ijj pe /uj ce ) 2 (vth,e/c) 2 = 0.02) they 
showed that only electrons are heated (accelerated), but 
not the ions. In our numerical simulations we see simi- 
lar behaviour: i.e. pre ferential acceleration of electrons 
(cf. Tsiklauri et al. 18]). In spite of this similarity, how- 
ever, one should make clear distinction between (a) phase 
mixing which essentially is a result of enhanced, due to 
plasma inhomogeneity, collisional viscous and ohmic dis- 
sipation and (b) resonant absorption which can operate 
in the collisionless regime as a result of generation of 
kinetic Alfven wave in the resonant layer, which subse- 
quently decays via Landau damping preferentially accel- 
erating electrons in the low (3 regime. 

(iii) They also showed that in the kinetic regime the 
total absorption rate is approximately the same as in 
the MHD. Our simulations seem to produce similar re- 
sults. We cannot comment quantitatively on this occa- 
sion, but at least the spatial from of the amplitude decay 
(oc exp(— Ax 3 )) is similar in both cases. 

Resuming aforesaid, we conjecture that the generated 
nearly parallel electrostatic fields found in our numerical 
simulations are due to the generation of kinetic Alfven 



waves that are created through interaction of initial shear 
Alfven waves with plasma density inhomogeneity, in a 
similar fashion as in resonant absorption described above. 
Further theoretical study is thus needed to provide solid 
theoretical basis for interpretation of our numerical sim- 
ulation results and to test the above conjecture. 
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FIG. 1: Contour (intensity) plots of phase mixed Alfven wave B z and E v components at instants: t = (15.62, 31.25, 54.69)/o; C i. 
Excitation source is at the left boundary. Because of periodic boundary conditions, left-propagating AW re-appears from the 
right side of the simulation box. Note how AW is progressively stretched because of differences in local Alfven speed. 
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FIG. 2: Left column: contour plots of generated nearly parallel to the external magnetic field electrostatic field E x at instants: 
t — (15.62, 31.25, 54. 69)/w c ;- Right column: ^-component of electron phase space at the same times. Note, that in order to 
reduce figure size, only electrons with V x > 0.15c were plotted. 
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FIG. 3: The distribution function of electrons at t = (dotted 
curve) and t — 54.69/aj C i (solid curve). 




FIG. 4: Two snapshots of J 
component at instants t = 
t = 46.87 /uJd (dotted line). 
0.056 exp [- (z/1250) 3 ]. 
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Dashed line represents fit 



